Shapefile Tally fo Choropleth

Overview

This document describes how to create a leaflet choropleth of the USA; a map of the USA split into regions (in this precise example, the 48 states of the contiguous USA) with each region shaded according to the count of data points contained within that region.

There are therefore two data sets required for this:

• Shape files: defining the choropleth regions
• Point data: sets of latitude and longtidues

The following libraries are required for this exercise:

library(maps)

## 
##  # maps v3.1: updated 'world': all lakes moved to separate new #
##  # 'lakes' database. Type '?world' or 'news(package="maps")'.  #

library(mapproj)
library(leaflet)
library(rgdal)

## Loading required package: sp

## rgdal: version: 1.1-10, (SVN revision 622)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 1.11.4, released 2016/01/25
##  Path to GDAL shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/gdal
##  Loaded PROJ.4 runtime: Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]
##  Path to PROJ.4 shared files: /Library/Frameworks/R.framework/Versions/3.3/Resources/library/rgdal/proj
##  Linking to sp version: 1.2-3

library(GISTools)

## Loading required package: maptools

## Checking rgeos availability: TRUE

## Loading required package: RColorBrewer

## Loading required package: MASS

## Loading required package: rgeos

## rgeos version: 0.3-19, (SVN revision 524)
##  GEOS runtime version: 3.4.2-CAPI-1.8.2 r3921 
##  Linking to sp version: 1.2-3 
##  Polygon checking: TRUE

## 
## Attaching package: 'GISTools'

## The following object is masked from 'package:maps':
## 
##     map.scale

library(sp)
library(plyr)

## 
## Attaching package: 'plyr'

## The following object is masked from 'package:maps':
## 
##     ozone

Load and Visualise Shape Files

Import the shapefiles using readOGR - obtained from https://www.census.gov/geo/maps-data/data/tiger-cart-boundary.html

## Load county borders
full_us_shapefiles <- readOGR(dsn = "data/shapefiles/", layer = "cb_2015_us_state_500k", verbose = F)
## Convert FIPS code to numeric
full_us_shapefiles$STATEFP <- as.character(full_us_shapefiles$STATEFP)
full_us_shapefiles$STATEFP <- as.numeric(full_us_shapefiles$STATEFP)

Regions in the USA are assigned FIPS codes, these may be used to identify those regions that belong to the contiguous USA using the data file “US-FIPS-Codes.csv”:

fips_codes <- read.csv("data/US-FIPS-Codes.csv", stringsAsFactors = F)
make_us_contiguous <- function(spatial_polgyon= NA){
  contiguous_fips_codes <- fips_codes[fips_codes$Contiguous.United.States. == "Y",]$STATE
  contiguous <- full_us_shapefiles[spatial_polgyon$STATEFP %in% contiguous_fips_codes,]
  # Drop unnecessary levels
  contiguous@data <- droplevels(contiguous@data)
  contiguous
}
contiguous_us_shapefiles <- make_us_contiguous(full_us_shapefiles)

Extract state names from FIPS file:

contiguous_fips_codes <- fips_codes[fips_codes$Contiguous.United.States. == "Y",]
contiguous_us_shapefiles$State_Name <- mapvalues(contiguous_us_shapefiles$STATEFP, 
                                                 from = contiguous_fips_codes$STATE, 
                                                 to = contiguous_fips_codes$STATE_NAME)

Visualise the shapefiles with leaflet:

map <- leaflet(data = contiguous_us_shapefiles) %>% addTiles()
map %>% addPolygons(
  stroke = FALSE,
  smoothFactor = 0.2,
  fillOpacity = 0.8,
  # fillColor = ~ pal(var),
  weight = 1
)

Count Points Per Location

Import the locations using read.csv and convert into a SpatialPointsDataFrame such that their intersection with the shapefiles may be easily calculated:

locations_df <- read.csv("data/lat-long-points.csv", stringsAsFactors = F)
locations_spdf <- SpatialPointsDataFrame(coords = locations_df,
                                         data = locations_df,
                                         proj4string = full_us_shapefiles@proj4string)

Use poly.counts to calculate the number of data points per polygon, noting that the polygons in the outputted tally are matched to the SpatialPolygonsDataframe object by row number - allowing the polygon tallies to be inserted into the SpatialPolygonsDataframe

## Example http://gis.stackexchange.com/a/113009
## Apply to my data!
contiguous_counts <- poly.counts(pts = locations_spdf, polys = contiguous_us_shapefiles)
contiguous_counts_df <- data.frame(contiguous_counts)
contiguous_us_shapefiles@data$Count.of.locations <- contiguous_counts_df$contiguous_counts

Create a colour palette using colorBin

palette <- colorBin(c("#cccccc",brewer.pal(5, "YlGnBu")), 
                    bins = c(0,1,5,10,20,50,350),
                    pretty = FALSE,
                    # na.color = "#cccccc",
                    alpha = TRUE
                    )

Generate choropleth:

region_labeller <- function(state_name = NA, number_of_points = NA){
  paste0(
    "<p>", state_name, "</p>",
    "<p>", number_of_points, "</p>"
  )
}

map <- leaflet(data = contiguous_us_shapefiles) %>% addTiles()
map <- map %>% addPolygons(
  stroke = TRUE,
  color = "#FFFFFF",
  smoothFactor = 0.2,
  fillOpacity = 0.8,
  fillColor = ~ palette(Count.of.locations),
  weight = 1,
  popup = ~region_labeller(state_name = State_Name, number_of_points = Count.of.locations)
)

map %>% addLegend(position = 'topleft', ## choose bottomleft, bottomright, topleft or topright
          colors = c("#cccccc",brewer.pal(5, "YlGnBu")),
          labels = c("0","1-5","5-10","10-20","20-50","50-350"),  ## legend labels (only min and max)
          opacity = 0.6,      ##transparency again
          title = "relative<br>amount")   ## title of the legend